Applicability of the hydrodynamic description of classical fluids 
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We investigate using numerical simulations the domain of applicability of the hydrodynamic de- 
scription of classical fluids at and near equilibrium. We find this to be independent of the degree 
of many-body correlations in the system; the range r c of the microscopic interactions completely 
determines the maximum wavenumber k max at which the hydrodynamic description is applicable 
by k max r c ~ 0.43. For the important special case of the Coulomb potential of infinite range, we 
show that the ordinary hydrodynamic description is never valid. 
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The equations of hydrodynamics [![, such as the 
Navier-Stokes equations, are certainly the most widely 
used framework for investigating the dynamics of fluids, 
including gases , liquids Q , plasmas [|| and nuclear 
matter [5(. Despite widespread use and successes, a num- 
ber of outstanding questions of both fundamental and 
practical importance remain regarding the conditions un- 
der which the hydrodynamic description holds. 

One ordinarily thinks of the hydrodynamic picture as 
applying only for wavenumbers k such that fclf < 1 
with If the mean free path and frequencies lo such that 
lo/u c <C 1 with uj c the mean collision frequency. These 
conditions, derived and already rather qualitative for a 
system governed by uncorrelated binary collisions (e.g. a 
dilute gas), become even more indeterminate when many- 
body correlations are present (as is the case, for exam- 
ple, in a liquid) because the concepts of mean free path 
and mean collision time cease to have a clear physical 
meaning. Thus the applicability of the hydrodynamic 
description certainly depends strongly on the thermody- 
namic conditions - e.g. the density n and temperature 
T - as well as the strength and range of the particles' 
mutual interactions. For instance, one expects that the 
description never applies on lengthscales smaller than the 
range r c of the potential (i.e. kr c > 1) - in other words, 
that the domain of validity will shrink as the range in- 
creases (such long range potentials are of particular im- 
portance in plasma physics). In fact, in the extreme case 
of r c = oo, the very existence of a hydrodynamic de- 
scription is a known but unsolved problem [6J. As well 
as leading to a deeper understanding of the emergence 
of macroscopic behaviour in interacting many-body sys- 
tems, it is also of significant practical importance to know 
exactly when hydrodynamics can be used to describe the 
behaviour of fluids, e.g. for analysing light and neutron 
scattering experiments [1, 0, Hi • 

In this Letter, we address with numerical simula- 
tions the question of the domain of applicability of the 
hydrodynamic description for fluids at or near equi- 
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librium as the level of many-body correlations in the 
system is varied. To this end, we consider a one- 
component system with a Yukawa interaction potential 
v(r) — q 2 exp(— r /r c )/r, where q 2 is the strength and 
r c the range of the potential. Although not possessing 
the short range features of conventional pair potentials 
used to describe normal liquids, the Yukawa potential 
is certainly suitable for investigating the long length- 
scale dynamics that concern the hydrodynamic descrip- 
tion. Additionally, it is commonly used in describing the 
screened interactions in plasmas [1, |T(| ■ What is more, 
since for r c — oo one recovers the Coulomb potential, 
we are able to use this model to answer the question of 
the existence of the hydrodynamic limit referred to pre- 
viously. Most importantly however, this model is very 
convenient here because it is known to be fully charac- 
terised by two dimensionless parameters only H - the 
reduced range r* = r c /a, where a = (47rn/3) is the 
average inter-particle spacing, and the coupling strength 
r = q 2 /(akbT), which itself characterises completely the 
degree of many-body correlations present in the system 
for a given range [11] . For a wide range of T and r* 
values, thus spanning states ranging from dilute gases to 
dense liquids [l2| . with short or long range microscopic 
interactions, we determine the length and time scales at 
which the hydrodynamic description breaks down. 

To accomplish this, we have computed with Molec- 
ular Dynamics (MD) simulations the dynamical struc- 
ture factor, S(k,uj), that is the Fourier transform in 
space and time of the density autocorrelation function, 
for a wide range of T (0.1,1,5,10,50,120,175) and r* 
(0.5,l,1.4,2,3.3,10,oo) values. S(k,w) contains complete 
information of the system dynamics at and near thermal 
equilibrium through the Fluctuation-Dissipation Theo- 
rem and is routinely measured in inelastic li ght and neu- 
tron scattering experiments (e.g. SB dill)- Three 
main difficulties are involved with MD calculation of 
S(k,oj). Firstly, for long range potentials (large r*), it 
is essential to include the Ewald summation; we com- 
pute this for all our r* values using the Particle-Particlc- 
Particle-Mesh method (Isj . Secondly, obtaining accu- 
rate MD data for S(k,uj) requires averaging the results 
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of a large number of simulations to improve statistics. 
Thirdly, in order to investigate the long wavelength dy- 
namics that concern the hydrodynamic description, very 
large scale simulations (a large number of particles N) are 
needed - the minimum reduced wavevector, (fca) m j n , at 
which the system dynamics can be determined using MD 
is oc TV -1 / 3 . These computational demands have made a 
thorough study like ours impractical before now. In our 
computation of S(k,uj), we average the results of fully 
25 simulations, each of duration 819. 2a;" 1 (the plasma 

frequency co p = y // 3q 2 ma 3 is the natural timescale for 
our system, where m is the particle mass), with up to 
500, 000 particles. Our complete analysis will be detailed 
in a forthcoming publication [l6j . 

Firstly, we consider the case of finite range interaction 
potentials (r* < oo). In this case the MD data can be 
compared to the result obtained from the linearised hy- 
drodynamic (Navier-Stokes) equations [U, E3 

S H (k,uj) _ 7- 1 2D T k 2 
S(k) ~ 7 oj 2 + (D T k 2 ) 2 ( ' 

ak 2 ok 2 
+ 1/7 [{uj + c s k) 2 + (ak 2 ) 2 + (uS - c s k) 2 + [ak 2 ) 2 \ ' 

where the static structure factor S(k) in Eq. (TTJ) is also 
determined from our MD simulations. Eq. (TTJ consists of 
a central (Rayleigh) peak representing a diffusive thermal 
mode and two Brillouin peaks at uj = ±c s fc correspond- 
ing to propagating sound waves. As illustrated in the top 
panel of Fig. [TJ at the smallest k value accessible to our 
MD simulations we find that the MD S(k, lu) can always 
(i.e. for all T and r*) be very accurately fitted to Eq. (fTJ) 
, thus giving numerical values for the thermal diffusiv- 
ity Dt, sound attenuation coefficient c, adiabatic sound 
speed c s and ratio of specific heats 7 that appear in the 
hydrodynamic description. When obtained in this way, 
these parameters are found to be in very good agreement 
with previous equation of state and transport coefficient 
calculations for the Yukawa model [3]. In particular, we 
find that 7 w 1 - that is, the Rayleigh peak at lu = in 
Yukawa fluids is negligible. In all cases, however, we find 
two Brillouin peaks, at 10 = ±c s fc, representing a damped 
sound wave. Fig. [2] shows the position of the Brillouin 
peak obtained from our MD simulations. We see that 
as the interaction potential becomes more long ranged, 
it is necessary to look at increasingly long lengthscales 
(small ka ) for the hydrodynamic description to be ap- 
plicable. Clearly in all cases, at some k value, which we 
denote by k max , the position Lj(k) of the Brillouin peak 
as computed by MD diverges from the linear relation. 

Quantitatively, we define k max as the minimum k value 
for which Lu(k)/(c s k) > 1.01. Using this criterion, for all 
values of the coupling T, we find that k max r c ~ 0.43. 

The k max obtained from the peak position is found 
to also characterise well the departure of the height and 
width of the Brillouin peak from the predictions of the 
hydrodynamic description (Fig. [3]). Therefore, k ma x is 
the maximum wavevector at which the hydrodynamic de- 
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FIG. 1: (color online) A sample of our MD results for S(k,u) 
(red dots) against S H (k,Lo) in Eq. |T} (black line) and when 
a 'mean field' is added (blue dashed line - bottom panel only) . 
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FIG. 2: (color online) Brillouin peak position w(k)/uj p as ob- 
tained from MD (open symbols), along with the correspond- 
ing linear relations ui — c s k (solid lines). 



scription of Eq. (fTJ) is applicable. As shown in Fig. |31 
beyond k max the height of the Brillouin peak decreases 
more slowly, and its width increases more slowly, than 
predicted by Eq. (fT]). Clearly however, the hydrody- 
namic description is valid for a relatively large range of k 
values, well beyond k = 0. In real space, we find that the 
lengthscale 2ir/k max is for all T values greater than the 
short-range correlation length over which the pair corre- 
lation function g[r) exhibits peaks and troughs @. It 
is remarkable that k max does not depend on T; indeed, 
one would intuitively expect the domain of validity of Eq. 
([T]) to increase as the system becomes more 'collisional' 
(i.e. with increasing T). We also note that providing 
k < k maxi the hydrodynamic approximation of Eq. (fTJ) 
for S(k,oj) is extremely accurate for all u> where S(k,u>) 
is not negligibly small; in this range, the Brillouin peaks 
exhaust the frequency sum-rules (see the top panel of 
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FIG. 3: (color online) Height and width of Brillouin peak as 
computed from MD (open circles), and the predictions of Eq. 
([1]) (solid lines). 

Fig. [Q). 

Much detailed work has been carried out to extend 
from macroscopic to microscopic lengthscales the do- 
main in which ordinary hydrodynamics applies (e.g. 
@, [H, E3])- Interestingly, we find that simply by adding 
to the usual stress tensor the mean field term one can 
account very well for the position of the Brillouin peak. 
Microscopically, this additional term stems from the in- 
clusion of a self-consistent 'mean field' or 'Vlasov' term 
- usually neglected because one considers lengthscales 
longer than the range of the potential - in the appropri- 
ate kinetic equation. By including the mean field term in 
the macroscopic equations, one obtains for the Yukawa 
model a modified expression for the position of the Bril- 
louin peak 19] 

where K — c 2 . — oJ^rj.. We note that for systems with 
7 = 1, which is a good approximation for the V and r* 
values considered here, the addition of the mean field 
does not change the hydrodynamic description of the 
height or width of the Brillouin peak (see [lj| for de- 
tails). As shown in Fig. 21 Eq. © gives a remarkably 
good description of the Brillouin peak position, even up 
to ka w 2 in most cases (although as shown in Fig. Q] 
the height and width of the peak does not always com- 
pare well with MD simulations). Indeed, this dramatic 
improvement is somewhat unexpected, since dynamics at 
these large wavevectors are not usually thought to be well 
described by macroscopic approaches. 

For finite r*, the mean field only begins to play a role 
when kr c > 0.43, i.e. when the range of the potential is 
large compared to the lengthscale of the density varia- 
tions. Therefore one may expect that for r* — oo, when 
the interaction potential is Coulombic [2(| , the mean field 
is important at all lengthscales (in this case, our criterion 
kmaxfc — 0.43 gives k max = !). To be sure, the pecu- 




FIG. 4: (color online) Brillouin peak position as obtained 
from MD (open circles), and the prediction of Eq. © (solid 
lines). 



liarity of the Coulomb potential is very well known - in 
this case the longitudinal waves are not low frequency 
sound waves as for r* < oo but instead high frequency 
plasma waves (uj w uj p ), even at k = 0. The resulting 
'plasmon' peak in S(k,uj), the position of which is illus- 
trated in Fig. [5] (red triangles), is certainly not described 
by the low-frequency hydrodynamic equations that lead 
to Eq. ([!} - one can indeed wonder why hydrodynamics 
should describe plasma oscillations at all. The difficulty 
here is underlined by a kinetic theoretical derivation of 
the hydrodynamic equations @: when proceeding with 
the Chapman-Enskog expansion of the appropriate ki- 
netic equation, the mean field term is usually treated 
as a small perturbation since in the small-gradient re- 
gion of interest to hydrodynamics the kinetic equation 
is always dominated by the collision term. In this case, 
however, the mean field term cannot be considered as 
small since its straightforward small-gradient expansion 
diverges with the characteristic Coulomb divergence (see 
@ and references therein). Based on this analysis, Baus 
and Hansen Q argued that only when the collisionality 
dominates the mean field, which they predicted would 
occur at sufficiently high coupling strength T, could a 
hydrodynamic description be expected. In this case the 
hydrodynamic description is identical to Eq. (QJ [2l| but 

with c s k replaced with uj p {l + -^jt) [22j . This macro- 
scopic description is known not to work for small T [f|; 
exactly how large T has to be for it to be applicable was 
left as an open question until now. We show here that in 
fact the hydrodynamic description is not valid at any T. 

Baus and Hansen [6] based the arguments outlined 
above on an exact formula for S(k, w), derived using gen- 
eralised kinetic theory , which at small k is given by [2l| 

S B (k,Lj) _ 5fc2 

S(k) ~ (u + u p (l + % a)) 2 + (f &) 2 



(w-w„(l + £ a)) 2 + (f by ' 
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r 


a 


2 / 2 

c a /u p 


b 


2a 


1 


0.895 


0.304 


0.192 


1.825 


5 


0.088 


-0.034 


0.109 


0.333 


10 


-0.009 


-0.080 


0.078 


0.235 


50 


-0.056 


-0.112 


0.032 


0.212 


120 


-0.062 


-0.127 


0.021 


0.349 


175 


-0.065 


-0.129 


0.009 


0.550 



TABLE I: Comparison of the generalised coefficients obtained 
by fitting the MD S(k,uj) at our smallest k value to Eq. Q 
with the analogous coefficients that appear in the hydrody- 
namic description. c s is determined from the internal energy 
fit given in [^|, while a is obtained from [l2| . 

where a and b are generalised coefficients with k and uj 
dependence. They were able to show that only at large 
r would it be possible for these coefficients to equal their 
macroscopic counterparts (of Eq. ([1])), c^/u^ and 2a 
respectively [6]. We have estimated a and b by fitting 
S(k,ui) at the smallest k value accessible to our MD sim- 
ulations to Eq. ([3]) - this gives a very good fit. However, 
as shown in Table[H the values obtained for a and b do not 
agree at all with their macroscopic counterparts, even at 
our highest coupling strength of T = 175, which is close 



to the freezing point of the system Q . For example, the 
width of the plasmon peak b does not even follow the 
same trend predicted by the hydrodynamic scaling 2a at 
our higher T values. From this we conclude that the com- 
bination of mean field and collisional effects means that 
the hydrodynamic description a la Navier Stokes is not 
valid for a Coulomb system at any coupling strength T. 

In summary, for finite range potentials, r* < oo , we 
find that the hydrodynamic description is (i) always valid 
at sufficiently long lengthscales where 'sufficiently long' is 
determined by the range of the potential (k max r c ~ 0.43) 
(ii) extremely accurate at these long lengthscales for all 
lo where S(k,uj) is not negligibly small (iii) not enlarged 
in its applicability as the level of many-body correlations 
in the system (i.e. T) is increased and (iv) is significantly 
extended in applicability by including a 'mean field' term 
in the macroscopic equations. For a Coulomb system, 
r* = oo, although the macroscopic approach correctly 
predicts the plasmon peak at k = 0, for k > the per- 
sistence of both mean field and collisional effects causes 
the ordinary hydrodynamic approach to fail. 
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